%% 设计巴特沃斯滤波器
T=2;
wp = 3;
ws = 7;
Wp = (2/T)*tan(wp/2);   % 设置归一化通带和阻带截止频率
Ws = (2/T)*tan(ws/2);
Rp = 1;
Rs = 40;  % 阻带最小衰减为40dB
% 双线性变换法
[N,Wc] = buttord(Wp,Ws,Rp,Rs,'s');   % 计算模拟滤波器的阶和3dB截止频率
[B,A] = butter(N,Wc,'low','s');      % 计算模拟滤波器系统函数分子、分母多项式系数
B
A
W = 0:0.001:50;
Hs = freqs(B,A,W); % 计算幅频响应
[Bz,Az] = bilinear(B,A,fs);      % 用双线性变换法法将模拟滤波器转换成数字滤波器
Hz = freqz(Bz,Az,W);     % 返回频率响应
phi1 = angle(Hs);
phi2 = angle(Hz);
Hs = 20*log10(abs(Hs));      % 幅度，单位为dB
% 画出滤波器的幅频响应和相频响应曲线
figure(13);
plot(W,Hs);
title('模拟滤波器的幅频响应');
xlabel('w');
ylabel('|Hw|');
grid on;
figure(14);
plot(W,phi1);
title('模拟滤波器的相频响应');
xlabel('w');
ylabel('angle(Hw)');
grid on;

% 画出数字滤波器的幅频响应和相频响应曲线
figure(15);
plot(W/pi,abs(Hz));
title('数字滤波器的幅频响应');
xlabel('w');
ylabel('|Hw|');
grid on;
figure(16);
plot(W,phi2);
title('数字滤波器的相频响应');
xlabel('w');
ylabel('angle(Hw)');
grid on;